TK!2
=v
#1
:n
x
:v
0
:s
i
:c
independent variable
:a
x
#2
:n
y1
:v
0
:s
i
:c
function value (numerical solution)
:a
y1
#3
:n
y2
:v
0
:s
i
:c
   "       "        "        "
:a
y2
#4
:n
yy1
:s
b
:c
function value (closed form)
:a
yy1
#5
:n
yy2
:s
b
:c
   "       "      "     "
:a
yy2
#6
:n
d1
:s
b
:c
difference
:a
d1
#7
:n
d2
:s
b
:c
    "
:a
d2
=u
=l
#1
:n
x
:c
independent variable
:v
0
.05
.1
.15
.2
.25
.3
.35
.4
.45
.4999999999999999
.6
.7
.7999999999999999
.8999999999999999
.95
1
E
#2
:n
Y
:c
names of unknown functions
:v
'y1
'y2
'y3
'y4
'y5
'y6
'y7
'y8
E
#3
:n
K
:c
names of Runge-Kutta coefficients
:v
'k1
'k2
'k3
'k4
'k5
'k6
'k7
'k8
E
#4
:n
k1
:v
13.31230446285602
13.98904825340202
14.02135526491138
14.73405756707261
E
#5
:n
k2
:v
13.60906626027465
14.24500530367558
14.28292296791739
14.95888824616458
E
#6
:n
yi
:v
6.696619749257822
6.642663058818072
E
#7
:n
ye
:v
7.397687512503392
7.356809207213942
E
#8
:n
y1
:c
numeric solution
:v
3
2.62434013671875
2.375340257024961
2.226371852490308
2.157610331867295
2.15444153522842
2.206255793213557
2.305536827251288
2.447175123300559
2.627952383431651
2.846156556780091
3.393934459288282
4.097801684569763
4.97760284366521
6.063795074150396
6.696619749257822
7.397179491479118
E
#9
:n
y2
:c
   "       "
:v
-7
-4.971506380208334
-3.39434817384637
-2.156194924648541
-1.171320127065752
-.374162947929487
.2855666376054379
.8466108239787327
1.33899732325934
1.78619753517375
2.206772493555465
3.024796560280331
3.884685974746178
4.854563985028177
5.992760607722579
6.642663058818072
7.356194817564949
E
#10
:n
yi0
:v
6.063795074150396
5.992760607722579
E
#11
:n
yh
:v
7.397169120864885
7.356158537082
E
#12
:n
yh1
:v
6.696619749257822
6.642663058818072
E
#13
:n
yh2
:v
7.397179491479118
7.356194817564949
E
#14
:n
DD
:v
.000001401968715926405
.000004931963311089608
E
#15
:n
yy1
:c
closed-form solution
:v
3
2.624315164525584
2.375302378921143
2.226328792505901
2.157566865037429
2.154400462309621
2.206218617632017
2.305504221182932
2.447147245217135
2.62792909167157
2.84613755087246
3.393883257539027
4.097759439721429
4.977587104201251
6.063814282271049
6.69665516219909
7.397229641807579
E
#16
:n
yy2
:c
     "         "
:v
-7
-4.9714060677241
-3.394195724883724
-2.15602113214359
-1.171143971943366
-.3739954957378435
.2857195314244758
.8467466526206573
1.339115661593796
1.786299189098466
2.206858938805384
3.025051583526627
3.884962075337657
4.854813705170566
5.992980192980528
6.642851562599982
7.356361927422938
E
#17
:n
d1
:c
differences  yy1 - y1
:v
0
-.00002497219316577315
-.00003787810381794543
-.00004305998440612413
-.00004346682986611228
-.0000410729187985126
-.0000371755815398167
-.00003260606835642221
-.00002787808342397113
-.00002329176008109002
-.00001900590763126431
-.00005120174925510312
-.00004224484833414976
-.00001573946395883041
.00001920812065314692
.00003541294126829797
.00005015032846067413
E
#18
:n
d2
:c
     "       yy2 - y2
:v
0
.000100312484233811
.0001524489626465453
.0001737925049507538
.0001761551223862234
.0001674521916434912
.0001528938190379359
.0001358286419246291
.0001183383344558475
.0001016539247160608
.00008644524991918701
.0002550232462956714
.0002761005914790893
.0002497201423894424
.000219585257948296
.0001885037819100432
.0001671098579887342
E
=f
#1
:n
RK4sa
:c
Fourth Order Runge-Kutta with selfadjusting step size
:t
Procedure
:i
ne,xz,eps
#2
:n
RK4step
:c
Fourth Order Runge-Kutta integration step
:t
Procedure
:i
ne,xi,h,y0,y
#3
:n
EQ
:c
1st order equations
:t
Procedure
:i
x,j
#1
:n
RK4sa
:z
" Notation:  ne   number of 1st order equations
"            xz   upper limit of the independent variable interval
"            eps  precision parameter
"            EQ   name of function with 1st order equations (see RK4step)
"            x    independent variable (list)
"            Y    list with names of lists representing function 
"                 values  y1, y2, ...
"            K    list with names of lists of RK coefficients  k1, k2, ...
"                 (see function subsheet for RK4step)

" Supporting function:  RK4step

" Description: This procedure implements a classical 4th order Runge-Kutta
" procedure for integration of ordinary differential equations.  The set of
" the 1st order differential equations must be available in a procedure
" function  EQ .  The initial values of x, y1, y2, ... must be set as the
" first element of the appropriate lists.  The starting step size is 1/10
" of the independent variable interval.  The steps are automatically 
" halved or doubled as required by the given precision parameter.
" Note: Reduction of the precision parameter below certain point may worsen
" rather than improve the solution accuracy because of accumulation of
" round-off errors from increased number of steps.

for j=1 to ne             " Initialization
   Yj:= 'Y[j]
   temp:= Yj[1]
   call blank(Yj)
   Yj[1]:= temp
next j
xa:= 'x[1]
call blank('x)
'x[1]:= xa
i2:= 1
h:= (xz-xa)/10             " starting with 10 steps
Loop:
   xi:= 'x[i2]
   for e=1 to ne
      'yi0[e]:= 'Y[e][i2]
   next e
   i1:= i2 + 1
   i2:= i1 + 1
   call RK4step(ne,xi,2*h,'yi0,'yh)    " taking a double-step
A: call listcopy ('yi0,'yi)
   call RK4step(ne,xi,h,'yi,'yh1)      " taking the 1st single step
   xi:= xi + h
   'x[i1]:= xi
   call listcopy('yh1,'yi)
   call RK4step(ne,xi,h,'yh1,'yh2)     " taking the 2nd single step
   for e=1 to ne
      DD:= 'yh[e] - 'yh2[e]
                                    " evaluating absolute or relative errors
      if abs('yh2[e])<eps then 'DD[e]:= abs(DD) else 'DD[e]:= abs(DD/'yh2[e])
      'Y[e][i1]:= 'yh1[e]
      'Y[e][i2]:= 'yh2[e]
   next e
   maxDD:= max('DD)
   if maxDD <= eps goto B
   call listcopy('yh1,'yh)
   xi:= xi - h
   h:= h/2                           " dividing current step by half
   goto A
B: 'x[i2]:= xi + h
   if (eps/maxDD) > 32 then h:= 2*h    " doubling size of next step
   if abs(xz-xi) > 2.5*abs(h) goto Loop
   h:= (xz-'x[i2])/2
   if abs(h) > max(1e-12,abs(xz)*1e-12) goto Loop
   'x[i2]:= xz
/E
E
#2
:n
RK4step
:z
" Description: This procedure performs an integration step in the Runge-
" Kutta algorithm as directed from the procedure RK4sa.

call listcopy(y0,'ye)
call EQ(xi,1)
for j=1 to 3
   H:= h/2
   if j=3 then H:=h
   x:=xi+H
   for e=1 to ne
      'ye[e]:= y0[e] + H*'K[e][j]
   next e
   call EQ(x,j+1)       " getting current values of RHS 
next j                  " of 1st order equations
for e=1 to ne
   y[e]:= y0[e] + dot('K[e],1,2,2,1)*h/6
next e
/E
E
#3
:n
EQ
:z
'K[1][j]:= .5*'ye[1] + 1.5*'ye[2]
'K[2][j]:= 6*'ye[1] - 4*'ye[2]
/E
E
=r
#1
:r
" ***  Solving ODE by Fourth Order Runge-Kutta method  ***
:s
C
#2
:r
" ***********  with selfadjusting step size  *************
:s
C
#4
:r
call RK4sa(2,1,.01)
#6
:r
yy1 = exp(2*x) + 2*exp(-5.5*x)      " closed-form solution
:s
C
#7
:r
yy2 = exp(2*x) - 8*exp(-5.5*x)           "     "     "
:s
C
#8
:r
d1 = yy1 - y1
:s
C
#9
:r
d2 = yy2 - y2
:s
C
#11
:r
" Notes:
:s
C
#13
:r
" This model is an example of using the Fourth Order Runge-Kutta method
:s
C
#14
:r
" for solving systems of ordinary differential equations.  The procedure
:s
C
#15
:r
" function RK4sa adjusts automatically the size ot the integration steps
:s
C
#16
:r
" as to maintain a predefined accuracy.
:s
C
#18
:r
" The given system of differential equations,
:s
C
#19
:r
"
:s
C
#20
:r
"         y1' =  .5*y1 + 1.5*y2   ;     y2' =  6*y1 - 4*y2
:s
C
#21
:r
"
:s
C
#22
:r
" is set in the procedure function  EQ  in terms of current values of
:s
C
#23
:r
" independent variable  x  and unknown functions  ye[1], ye[2] .
:s
C
#25
:r
" The values of independent variable and unknown functions are generated
:s
C
#26
:r
" during the solution in lists  x  and  y1, y2  respectively.  The initial
:s
C
#27
:r
" values of  x, y1, y2  have to be typed in as the first elements of these
:s
C
#28
:r
" lists, most conveniently in the interactive table  Results .  Once this is
:s
C
#29
:r
" done, a call to function  RK4sa  triggers the integration and filling the
:s
C
#30
:r
" lists with the remaining values.  It is fun to watch this happening when
:s
C
#31
:r
" in the interactive table with 'Display Intermediate Values' on 'Yes'
:s
C
#32
:r
" (although switching it off speeds up the solution process some 10 times).
:s
C
#33
:r
" It is especially interesting to watch how the size of the integration
:s
C
#34
:r
" step (i.e. spacing in list x) changes as the solution proceeds.  The
:s
C
#35
:r
" graphs of the resulting functions may be displayed by pressing F7.
:s
C
#37
:r
" The linear differential equations in this example have a simple closed-
:s
C
#38
:r
" form solution (see 2nd and 3rd rule above) so that the accuracy of the
:s
C
#39
:r
" numerical solution may be checked and the effect of changing the precision
:s
C
#40
:r
" parameter (the last argument in RK4as) may be investigated.  In order to
:s
C
#41
:r
" do so the 1st rule must be cancelled, all other rules uncancelled, and 
:s
C
#42
:r
" list solving [F10] applied.
:s
C
=p
#1
:n
Results
:t
Line chart
:d
EGA
:i
""
:s
Yes
:z
X-axis
:g
No
:l
None
:h
""
:v
""
:x
x
:u
#1
:y
y1
:f
1
:p
Lines
:c
"*"
#2
:y
y2
:f
1
:p
Lines
:c
"*"
E
=t
#1
:n
Results
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
""
:u
#1
:l
x
:w
5
:h
""
#2
:l
y1
:w
10
:h
""
#3
:l
y2
:w
10
:h
""
#4
:l
yy1
:w
10
:h
""
#5
:l
yy2
:w
10
:h
""
#6
:l
d1
:n
dform
:w
10
:h
""
#7
:l
d2
:n
dform
:w
10
:h
""
E
=n
#1
:n
dform
:d
15
:f
3
:^
Scientific
:e
"."
:m
""
:p
Zero
:j
Right
:s
- Only
:q
0
:x
0
:l
""
:t
""
:v
"0.000E+0"
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tr,f,0,1,0,1,0
